      subroutine fence_box(ncell, ijkloc, iphead, x_phys, y_phys, z_phys, up, &
         vp, wp, xl, xr, yb, yt, ze, zf, wl, wr, wb, wt, we, wf, killer, npart) 
! 	A routine to impose BC on the particles
! 	1= reflection
!	2= periodic
!	other = removal
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ncell 
      integer , intent(in) :: wl 
      integer  :: wr 
      integer , intent(in) :: wb 
      integer  :: wt 
      integer , intent(in) :: we 
      integer  :: wf 
      integer , intent(in) :: npart 
      real(double) , intent(in) :: xl 
      real(double) , intent(in) :: xr 
      real(double) , intent(in) :: yb 
      real(double) , intent(in) :: yt 
      real(double) , intent(in) :: ze 
      real(double) , intent(in) :: zf 
      integer , intent(in) :: ijkloc(*) 
      integer , intent(in) :: iphead(*) 
      integer , intent(out) :: killer(*) 
      real(double) , intent(inout) :: x_phys(0:npart) 
      real(double) , intent(inout) :: y_phys(0:npart) 
      real(double) , intent(inout) :: z_phys(0:npart) 
      real(double) , intent(inout) :: up(0:npart) 
      real(double) , intent(inout) :: vp(0:npart) 
      real(double) , intent(inout) :: wp(0:npart) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk, np, m, modd 
      real(double) :: zm2, xx, x1, fact, yy, y1, zz, z1 
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!-----------------------------------------------
      integer , external :: cvmgp 
      real(double) , external :: ccvmgz 
!-----------------------------------------------
!
!
!	1 - riflessione
!	2 - periodicita
!	other - rimozione
!
!
      zm2 = ze + zf 
!
      killer(:ncell) = 0 
!
!     x direction
!
      xx = xr - xl 
 
      if (wl == 1) then 
 
!dir$ ivdep
         do n = 1, ncell 
            ijk = ijkloc(n) 
            np = iphead(ijk) 
 
            x1 = x_phys(np) - xl 
            m = ifix(abs(x1)/xx) 
            m = m*cvmgp(1,-1,x1) + cvmgp(0,-1,x1) 
            x_phys(np) = x_phys(np) - m*xx 
            modd = abs(mod(m,2)) 
            x_phys(np) = x_phys(np)*(1 - 2*modd) + modd*xx 
!
!     implements Pritchett's bc
!
            fact = ccvmgz(0.,1.,m) 
!
!     diagnostics
!
!     rrr=x1/xx
!     if(fact.eq.0.) then
!     if(rrr.ge.0.0.and.rrr.le.1.0) then
!c    ok
!     else
!     write(*,*)'ERRRORE ERRORE',abs(x1)/xx,fact
!     write(60,*)'ERRRORE ERRORE',abs(x1)/xx,fact
!     end if
!     else
!     if(rrr.lt.0.0.or.rrr.gt.1.0) then
!c    ok
!     else
!     write(*,*)'ERRRORE ERRORE',abs(x1)/xx,fact
!     write(60,*)'ERRRORE ERRORE',abs(x1)/xx,fact
!     end if
!     end if
!
!     if(fact.eq.1.)write(*,*)'pre',up(np),wp(np),z_phys(np),fact
            up(np) = (1. - fact)*up(np) - up(np)*fact 
!      uncomment the next two lines for Pritchett's bc
!     wp(np)=(1.-fact)*wp(np)-wp(np)*fact
!     z_phys(np)=(1.-fact)*z_phys(np)+(zm2-z_phys(np))*fact
!     if(fact.eq.1.)write(*,*)'post',up(np),wp(np),z_phys(np)
         end do 
 
      else if (wl == 2) then 
 
!dir$ ivdep
         do n = 1, ncell 
            ijk = ijkloc(n) 
            np = iphead(ijk) 
 
            x1 = x_phys(np) - xl 
            m = ifix(abs(x1)/xx) 
            m = m*cvmgp(1,-1,x1) + cvmgp(0,-1,x1) 
            x_phys(np) = x_phys(np) - m*xx 
 
         end do 
 
      else 
 
!dir$ ivdep
         do n = 1, ncell 
            ijk = ijkloc(n) 
            np = iphead(ijk) 
 
            x1 = x_phys(np) - xl 
            m = ifix(abs(x1)/xx) 
            m = m*cvmgp(1,-1,x1) + cvmgp(0,-1,x1) 
            x_phys(np) = max(min(x_phys(np),xr),xl) 
            killer(n) = killer(n)+ abs(m) 
 
         end do 
 
      endif 
 
!
!     y direction
!
      yy = yt - yb 
 
      if (wb == 1) then 
 
!dir$ ivdep
         do n = 1, ncell 
            ijk = ijkloc(n) 
            np = iphead(ijk) 
 
            y1 = y_phys(np) - yb 
            m = ifix(abs(y1)/yy) 
            m = m*cvmgp(1,-1,y1) + cvmgp(0,-1,y1) 
            y_phys(np) = y_phys(np) - m*yy 
            modd = abs(mod(m,2)) 
            y_phys(np) = y_phys(np)*(1 - 2*modd) + modd*yy 
 
            fact = ccvmgz(0.,1.,m) 
            vp(np) = (1. - fact)*vp(np) - vp(np)*fact 
 
         end do 
 
      else if (wb == 2) then 
 
!dir$ ivdep
         do n = 1, ncell 
            ijk = ijkloc(n) 
            np = iphead(ijk) 
 
            y1 = y_phys(np) - yb 
            m = ifix(abs(y1)/yy) 
            m = m*cvmgp(1,-1,y1) + cvmgp(0,-1,y1) 
            y_phys(np) = y_phys(np) - m*yy 
 
         end do 
 
      else 
 
!dir$ ivdep
         do n = 1, ncell 
            ijk = ijkloc(n) 
            np = iphead(ijk) 
 
            y1 = y_phys(np) - yb 
            m = ifix(abs(y1)/yy) 
            m = m*cvmgp(1,-1,y1) + cvmgp(0,-1,y1) 
            y_phys(np) = max(min(y_phys(np),yt),yb) 
            killer(n) = killer(n)+ abs(m) 
 
         end do 
 
      endif 
!
!     z direction
!
      zz = zf - ze 
 
      if (we == 1) then 
 
!dir$ ivdep
         do n = 1, ncell 
            ijk = ijkloc(n) 
            np = iphead(ijk) 
 
            z1 = z_phys(np) - ze 
            m = ifix(abs(z1)/zz) 
            m = m*cvmgp(1,-1,z1) + cvmgp(0,-1,z1) 
            z_phys(np) = z_phys(np) - m*zz 
            modd = abs(mod(m,2)) 
            z_phys(np) = z_phys(np)*(1 - 2*modd) + modd*zz 
 
            fact = ccvmgz(0.D0,1.D0,m) 
            wp(np) = (1. - fact)*wp(np) - wp(np)*fact 
 
 
         end do 
 
      else if (we == 2) then 
 
!dir$ ivdep
         do n = 1, ncell 
            ijk = ijkloc(n) 
            np = iphead(ijk) 
 
            z1 = z_phys(np) - ze 
            m = ifix(abs(z1)/zz) 
            m = m*cvmgp(1,-1,z1) + cvmgp(0,-1,z1) 
            z_phys(np) = z_phys(np) - m*zz 
 
         end do 
 
      else 
 
!dir$ ivdep
         do n = 1, ncell 
            ijk = ijkloc(n) 
            np = iphead(ijk) 
 
            z1 = z_phys(np) - ze 
            m = ifix(abs(z1)/zz) 
            m = m*cvmgp(1,-1,z1) + cvmgp(0,-1,z1) 
            z_phys(np) = max(min(z_phys(np),zf),ze) 
            killer(n) = killer(n)+ abs(m) 
 
         end do 
 
      endif 
 
 
      return  
      end subroutine fence_box 


 
      subroutine map_box(ncell, ijkloc, iphead, x_phys, y_phys, z_phys, x_nat, &
         y_nat, z_nat, ibar, jbar, kbar, dx, dy, dz, npart, nu_len, nu_comm) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ncell 
      integer , intent(in) :: ibar 
      integer , intent(in) :: jbar 
      integer , intent(in) :: kbar 
      integer , intent(in) :: npart 
      real(double) , intent(in) :: dx 
      real(double) , intent(in) :: dy 
      real(double) , intent(in) :: dz 
      real(double) , intent(in) :: nu_len 
      logical , intent(in) :: nu_comm 
      integer , intent(in) :: ijkloc(*) 
      integer , intent(in) :: iphead(*) 
      real(double) , intent(in) :: x_phys(0:npart) 
      real(double) , intent(in) :: y_phys(0:npart) 
      real(double) , intent(in) :: z_phys(0:npart) 
      real(double) , intent(inout) :: x_nat(0:npart) 
      real(double) , intent(inout) :: y_nat(0:npart) 
      real(double) , intent(inout) :: z_nat(0:npart) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: nn, ic, n, ijk, np 
      real(double) :: r, diz 
!-----------------------------------------------
!
      r = kbar*dz 
      nn = kbar + 1 
      ic = (nn - 1)/2 + 2 
      diz = nn/nu_len 
!
!dir$ ivdep
      if (.not.nu_comm) then 
         do n = 1, ncell 
            ijk = ijkloc(n) 
            np = iphead(ijk) 
            if (np <= 0) cycle  
            x_nat(np) = x_phys(np)/dx + 2. 
            y_nat(np) = y_phys(np)/dy + 2. 
!
!     uniform grid
!
            z_nat(np) = z_phys(np)/dz + 2. 
!
            x_nat(np) = max(2.,min(float(ibar + 2) - 1.E-10,x_nat(np))) 
            y_nat(np) = max(2.,min(float(jbar + 2) - 1.E-10,y_nat(np))) 
            z_nat(np) = max(2.,min(float(kbar + 2) - 1.E-10,z_nat(np))) 
 
         end do 
      else 
         do n = 1, ncell 
            ijk = ijkloc(n) 
            np = iphead(ijk) 
            if (np <= 0) cycle  
            x_nat(np) = x_phys(np)/dx + 2. 
            y_nat(np) = y_phys(np)/dy + 2. 
!
!     non-uniform grid in z
!
!     z_nat(np)=diz*asinh(z_phys(np)/R*(sinh((NN+1-ic)/diz)
!     &           -sinh((2-ic)/diz))+sinh((2-ic)/diz))+ic
            z_nat(np) = z_phys(np)/r*(sinh((nn + 1 - ic)/diz) - sinh((2 - ic)/&
               diz)) + sinh((2 - ic)/diz) 
            z_nat(np) = diz*log(z_nat(np)+sqrt(z_nat(np)**2+1.)) + ic 
!
            x_nat(np) = max(2.,min(float(ibar + 2) - 1.E-10,x_nat(np))) 
            y_nat(np) = max(2.,min(float(jbar + 2) - 1.E-10,y_nat(np))) 
            z_nat(np) = max(2.,min(float(kbar + 2) - 1.E-10,z_nat(np))) 
 
         end do 
      endif 
 
      return  
      end subroutine map_box 
